Project Goal: To compare the preliminary SWOT L2_HR_Raster product with a widely adopted surface water map (Global Surface Water Occurrence – GSWO)
Background: The Surface Water and Ocean Topgraphy (SWOT) satellite launched in December of 2022, and will vastly enrich scientific knowledge of Earth’s water cycle. The Ka-band Radar Interferometer (KaRIn) on board sends radar pulses Earth, and collects those pulses with two antennae. As the radar pulses return they reach each antenna with a different phase, because they traveled different distances (path lengths). Using the principles of physics and geometry underlying electromagnetic waves, NASA’s algorithms can measure the elevation and presence of surface water. SWOT will measure global inland surface water with high vertical (10cm) and horizontal (15cm) accuracy twice (on average) every 21 days.
Here’s a fantastic animation NASA made showing SWOT in action! You can see the orbital path, measurement footprint and radar pulses. (https://swot.jpl.nasa.gov/resources/142/swot-global-coverage/)
The vast majority of SWOT data is still not publicly available; however, NASA released a small preliminary beta version for scientists to familiarize themselves with the dataset’s limitations and opportunities. The beta release is also a chance for users to build workflows and learn about the complex file structures. In this analysis, I’ll explore a preliminary scene from SWOT and compare it to a widely adopted surface water map (GSWO).
Area of Study: I selected a SWOT image from Oregon, which contains the Columbia River and the area around Portland/Vancouver.
Read the SWOT in its ncdf format.
Read the Global Surface Water Occurance dataset (GSWO). Clip, project, and resample GSWO to match SWOT’s extent, projection and resolution.
Compare the SWOT and GSWO data sets.
Modify the SWOT data set to better match the GSWO data.
# The SWOT path names contain rich metadata such as observation time, cycle #, UTM zone, etc. After downloading, I changed the file name for simplicity.
swot_path <- paste0(data_raw_path, 'SWOT1.nc')
swot_nc <- nc_open(swot_path)
# Pulls the water
x <- ncvar_get(swot_nc, 'x')
y <- ncvar_get(swot_nc, 'y')
wtr_area_frac_array <- ncvar_get(swot_nc, 'water_area')
# This function gets ncdf file's meta data, which is handy in later steps.
swot_nc_atts <- ncatt_get(swot_nc, 0)min_x <- swot_nc_atts$geospatial_lon_min
min_y <- swot_nc_atts$geospatial_lat_min
max_x <- swot_nc_atts$geospatial_lon_max
max_y <-swot_nc_atts$geospatial_lat_max
bbox <- c(xmin = min_x, xmax = max_x, ymin = min_y, ymax = max_y)
# Create a spatial extent object from the bounding box
bbox_extent <- extent(bbox)
# Clip the raster using the bounding box
gsw_occurance_bbox <- crop(gsw_occurance, bbox_extent)
# Clean up memory
rm(gsw_occurance)# Let's clip the GSWO data to the SWOT image's footprint
# The SWOT data is projected to UTM10N; however, these bounds are in WSG:84 coordinates.
# Clipping the GSWO data before projecting, will reduce the computational load.
swot_lf_lon <- swot_nc_atts$left_first_longitude
swot_lf_lat <- swot_nc_atts$left_first_latitude
swot_ll_lon <- swot_nc_atts$left_last_longitude
swot_ll_lat <- swot_nc_atts$left_last_latitude
swot_rf_lon <- swot_nc_atts$right_first_longitude
swot_rf_lat <- swot_nc_atts$right_first_latitude
swot_rl_lon <- swot_nc_atts$right_last_longitude
swot_rl_lat <- swot_nc_atts$right_last_latitude
coords <- matrix(c(swot_lf_lon, swot_lf_lat, # left_first
swot_ll_lon, swot_ll_lat, # left_last
swot_rl_lon, swot_rl_lat, # right_last
swot_rf_lon, swot_rf_lat, # right_first
swot_lf_lon, swot_lf_lat), # closing the polygon by repeating the first point
byrow = TRUE, ncol = 2)
poly <- vect(coords, type = "polygons")
crs(poly) <- "EPSG:4326"
gsw_occurance_tiled <- mask(gsw_occurance_bbox, poly, inverse = FALSE)## Reproject the GSWO data to match SWOT’s projection. Resample the GSWO bringing it from 30m resolution to SWOT’s 100m resolution. ##
## Warning: [rast] CRS do not match
wtr_vector <- na.omit(as.vector(wtr_area_frac_array))
wtr_frac_df <- data.frame(value = wtr_vector)
# Quick histogram for the SWOT
pixel_histogram <- ggplot(data = wtr_frac_df,
mapping = aes(x = value)) +
geom_histogram(bins = 500) +
xlim(-10, 100000) +
geom_vline(xintercept = 10000, color = 'red', linewidth = 1, linetype = "dashed") +
theme_bw() +
labs(title = "Distribution of water area pixels",
y = "Count (log scale)",
x = "Water pixel value (square meters)") +
scale_y_log10()
(pixel_histogram)## Warning: Removed 17217 rows containing non-finite values (`stat_bin()`).
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 31 rows containing missing values (`geom_bar()`).
The mask pixel values shouldn’t exceed 10,000 square meters, and should
never be less than zero. According to SWOTs documentation, erroneous
pixels values are included for users to process at their discretion. The
red line on this plot marks the 10,000 square meters cutoff.